Weak plaquette valence bond order in the S = 1/2 honeycomb Ji — J2 Heisenberg model 
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Using the density matrix renormalization group, we investigate the 5=1/2 Heisenberg model 
on the honeycomb lattice with first- (Ji) and second-neighbor (J2) interactions. We are able to 
study long open cylinders with widths up to 12 lattice spacings. For J2/J1 near 0.3, we find an 
apparently paramagnetic phase, bordered by an antiferromagnetic phase for J2 < 0.26 and by a 
valence bond crystal for J2 > 0.36. The longest correlation length that we find in this intermediate 
phase is for plaquette valence bond (PVB) order. This correlation length grows strongly with 
cylinder circumference, indicating either quantum criticality or weak PVB order. 

PACS numbers: 



Progress in finding realistic model quantum Hamilto- 
nians with spin-liquid (SL) ground states has accelerated 
dramatically in the last two years, almost 40 years since 
Anderson first proposed a resonating valence bond (RVB) 
state as a possible ground state of the triangular Heisen- 
berg model [l| . One key recent advance was the discovery 
using the density matrix renormalization group (DMRG) 
of a gapped SL ground state in the spin- 1/2 kagome 
Heisenberg antiferromagnet 0, Spin liquid phases 
have been suggested for various other models, such as 
the half-filled honeycomb Fermi-Hubbard model [4j and 
the square lattice spin-1/2 Heisenberg antiferromagnet 
with second-neighbor (J2) interactions [5, 6]. However, 
some skepticism has been expressed about the evidence 
for spin liquids in the latter two models 0, Q ■ 

The main defining feature of a quantum spin liquid is 
the absence of any spontaneously broken symmetry, par- 
ticularly either magnetic or valence-bond order. Frus- 
tration, which discourages order, is a key ingredient of 
models potentially containing spin liquid phases. Spin 
liquids arise in several analytic treatments and exactly 
solvable, simplified, but less realistic models [9]. A key 
feature distinguishing types of spin liquids is the pres- 
ence or absence of a gap to all excitations. The kagome 
Heisenberg spin liquid is found to be gapped. To satisfy 
the Lieb-Schultz-Mattis theorem, gapped spin liquids for 
models with a net half-integer spin per unit cell must have 
"hidden" topological degeneracies in the thermodynamic 
limit, which depend on the topology of the system. The 
simplest possibility is a Z2 spin liquid. Since local mea- 
surements cannot identify Z2 or other topological orders, 
it is challenging to identify its presence in a numerical 
study. The degeneracies characteristic of a 2D gapped 
Z2 spin liquid have not been accessible for the system 
sizes studied to date. Odd-width cylinders spontaneously 
dimerize in a pattern that is characteristic of a quasi-one- 
dimensional system 0, Q . Another key feature of a Z2 
spin liquid is the presence of a — In 2 constant term cor- 
rection to the linear growth of the entanglement entropy 
with subsystem perimeter. This term has now been mea- 



sured in the nearest-neighbor kagome system [3| and also 
in the kagome system with next-nearest-neighbor interac- 
tion J2 pJI], where for J2 =0.1 the gaps are large and the 
entanglement entropy correction term can be measured 
particularly precisely. Thus, there is now solid evidence 
that the ground state of the kagome spin-1/2 antiferro- 
magnet is a gapped Z2 spin liquid. 

In this paper, we examine the spin-1/2 Heisenberg an- 
tiferromagnet on the honeycomb lattice (see Fig. [TJa)) 
with Hamiltonian 



H — Ji ^ Si ■ Sj + J2 ^ Si ■ Sj 



(1) 



where the sum over runs over nearest-neighbor 

pairs of sites and the sum ((«, j)) runs over next-nearest- 
neighbors. We take Ji = 1 (antiferromagnetic (AF)) 
and consider only J2 > 0. Our work follows other stud- 
ies of this and similar Heisenberg models, motivated by 
the Hubbard model results 
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. Most of these stud- 
ies report a nonmagnetic phase near J2/J1 ^ 0.2 — 0.4, 
and we agree, but the results are in general disagreement 
on the range and nature of this phase. A variational 
Monte Carlo study indicated a spin liquid in the range 
0.08 to 0.3 [20]. A combination of exact diagonalization 
and valence-bond treatment reported a plaquette valence 
bond (PVB, see Fig. [Hb)) crystal in the ran ge .2-0.4 
11 1 . Exact diagonalization on small lattices [1^, and 
the coupled-cluster method both suggest plaquette va- 
lence bond (PVB) order 13]. Earlier work reported that 
dimer correlations aren't strong enough for PVB order 
[l6| . Functional renormalization group work claims this 
phase has weak dimer and plaquette response 15| . Vari- 
ational entangled plaquette states suggest that none of 
the order parameters remain nonzero [l9[. Other theo- 
retical work has focused on the possibility of a Z2 SL on 
the honeycomb lattice and on phase transitions between 
Neel and staggered valence-bond crystal (SVBC) phases 
[21-23]. 

Here we report that the ground state displays an ap- 
parently paramagnetic phase for 0.26 ^ J2 S 0.36. For 
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FIG. 1: (color online) (a) One hexagon of the honeycomb lattice showing Ji and J2 interactions. The six sides of the hexagon 
are Ji bonds and the six J2 bonds are shown as (blue) dashed lines, (b) Illustration of the pattern of plaquette valence bond 
(PVB) order. The spin-spin correlations differ between the bonds shown thin (black) vs. thick (red). (c)Three phases on 
the YC6-0 cylinder with 300 sites and a gradient of J2. The widths of lines are proportional to \{Si ■ Sj)\. We also show the 
second-neighbor spin correlation, but only when its value is negative. The arrows represents the values of {Sf) at each site. 
The scales for correlations and magnetizations are indicated. We remove several bonds in one zigzag row of the cylinder, so 
that the AF spin pattern is more clearly seen. The two vertical lines at J2 — 0.26 and 0.36 indicate estimates of the phase 
boundaries. 



J2 < 0.26, we find the usual two-sublattice AF phase. 
For J2 ^ 0.36, we find a staggered valence bond crys- 
tal (SVBC). We have studied in some detail the system 
with J2 = 0.3, deep within the intermediate phase that 
is neither AF nor SVBC. We examine various correla- 
tion functions of this ground state on various cylinders. 
The longest correlation length that we find is for PVB 
order. The strong growth of the PVB correlation length 
with cylinder circumference indicates that the system is 
either near a quantum critical point or may have weak 
long-range PVB order. 

We use cylindrical (C) boundary conditions with open 
ends for our DMRG [Ij, IHj calculations. We label the 
cylinders either XCM-N or YCM-N. The labels X or Y 
indicated whether a first-neighbor bond is oriented hor- 
izontally (X) or vertically (Y). For XC cylinders, M is 
the number of sites along a zigzag vertical column and 
N means that the periodic boundary conditions are con- 
nected with a shift of N zigzag columns to the left or right. 
For YC cylinders, M is the number of zigzag horizontal 
rows and N means the connection has a shift by N sites 
along a zigzag row. For example, in the XC8-0 cylinder, 
one set of edges of each hexagon lie along the X direction, 
and there are 8 sites along the zigzag columns, which are 
connected periodically. So the circumference is C = 4-\/3 
lattice spacings. For the YC4-0 cylinder, one set of edges 
of each hexagon lie along the Y direction and the cylin- 
der is connected periodically along the Y direction with 
circumference 6 lattice spacings. For the XC9-1 cylin- 
der, the connection has a horizontal shift of one zigzag 
column, producing a circumference of C = 3^/7. 

In Fig. [TJc), we present the ground state of a single 
system which gives an overview of the entire phase dia- 
gram. For this long YC6-0 cylinder, J2 is uniform along 
the vertical direction, but varies linearly with the hori- 



zontal position from J2 = at the left edge to J2 = 0.5 at 
the right edge. To make the AF order visible, we apply 
a staggered field at the left end of the cylinder. As J2 
increases along the cylinder, the AF order decreases, be- 
coming negligible in the intermediate phase. We will dis- 
cuss this intermediate state in detail below. The SVBC 
phase appears clearly for J2 > 0.36. This SVBC phase 
has strong first-neighbor correlations along the vertical 
direction with strong horizontal second-neighbor correla- 
tions connecting them to form "ladders" . Below we will 
determine the phase boundaries of the AF and SVBC 
phases more accurately. 

First we estimate the boundary of the AF phase. One 
technique to determine magnetic order parameters using 
DMRG is to put strong ordering fields on the edges of 
an open cylinder, and adjust the aspect ratio Ly/L^ to 
to minimize the finite size effects [26|. For both square 
and triangular spin- 1/2 Heisenberg antiferromagnets, an 
aspect ratio near 1.7 ~ 1.9 is found to minimize the finite 
size effects. For XCM-0 cylinders with M columns, the 
aspect ratio is a/3, which we use. For J2 = 0, we deter- 
mine that (Sz) = 0.27, which is close to the value de- 
termined using Monte Carlo in the thermodynamic limit 
(Sz) = 0.2677(6) [23|. With J2 increasing, we find that 
the magnetization reduces to near zero for J2 = 0.26 in 
Fig. [5] The various cluster sizes all point to the phase 
transition near 0.26. This phase transition point is larger 
than the classical limit value of J2 = |. Ref. [2^ claims 
the Neel order disappears at 0.08, however the value we 
find here is more consistent with other studies which give 
J2 = 0.2 [ni[l4|. 



To estimate the boundary of the SVBC phase we study 
several XC cylinders with J2 varying from 0.28 to 0.40 
using the method in Fig. [T] These results show the 
SVBC phase for J2 > 0.36. We also use the entangle- 
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FIG. 2: (color online) The staggered magnetization at the 
center of the cylinder versus J2 for various XC cylinders. The 
inset shows how the local magnetization decays from the edge 
of a long XCIO cylinder for various values of J2. 
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FIG. 3: (color online) The entanglement entropy and its 
derivative versus J2 for different XC cylinders near the tran- 
sition in to the SVBC phase. 



ment entropy and its first derivative virith respect to J2 
to estimate the phase transition point. We make vertical 
cuts between zigzag columns, the dividing line between 
the two parts of the system bisecting a column of hor- 
izontal bonds, and measure the entanglement entropy. 
As seen in Fig. |31 the entropy drops in going from the 
intermediate phase to the SVBC. The derivative of the 
entropy shows a peak around 0.37 for the XC8-0 cylin- 
der, and around 0.36 for the wider XClO-0 and XC12-0 
cylinders. In addition, the height of this peak increases 
with the system width, as expected for a peak indicating 
a phase transition. |28l430| 

On XC cylinders, we find the SVBC state has strong 
first- and second-neighbor correlations along diagonal di- 
rections, forming diagonally oriented ladders. Thus there 
are two degenerate diagonal SVBC states on an XC 
cylinder, whereas for YC cylinders there is only the one 
vertical SVBC pattern (Fig. 2). For an infinite two- 
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FIG. 4; (color onhne) Weak PVB order at J2 = 0.3 on (a) 
YC7-3, (b) XC9-1 and (c) XC12-0 cylinders. The widths of 
the lines are proportional to \{Si ■ Sj) -|-0.32| for all the plots, 
with blue solid(red dashed) lines for negative(positive) val- 
ues. The bond strength scale is indicated at the top. The 
three different sublattices of plaquettes are labeled as A, B, 
C. In these figures we keep m = 6000 states in our DMRG 
calculation. The truncation error for XC9-1 is smaller than 
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while it is near 10~^ for the wider YC7-3 and XC12-0 



cylinders. 



dimensional system all 3 of these SVBC ground states 
would be degenerate by rotational symmetry. In the clas- 
sical limit, for large J2 values, the ground state is a spin 
spiral state. However, quantum fiuctuations are strong 
enough to melt the spiral order and form the SVBC, in 
agreement with rcf. [17]. 

In the rest of this paper, we focus specifically on 
J2 — 0.3 inside the intermediate phase. To measure mag- 
netic correlations, we apply "pinning" magnetic fields at 
one end of the cylinder, and measure the resulting mag- 
netization pattern. Unlike in the AF phase, the induced 
magnetization decays exponentially from the end of the 
cylinder with a decay length of 2 to 3 lattice spacings for 
various cylinders. We further check the response to a lo- 
cal magnetic field applied to a spin at the cylinder center. 
This local magnetic field response is quite short ranged. 
It only infiuences its nearby surrounding sites, as opposed 
to generating a large region of staggered magnetization 
in the AF phase. 

Refs. 



Ill , |12| suggest that the intermediate phase is a 



PVB phase with long range dimer-dimer correlations. To 
investigate PVB ordering we study cylinders with peri- 
odic boundary conditions that are compatible with PVB 
order, including YC4-0, YC6-0, YC7-3, YC8-0, XC6-0, 
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FIG. 5: (color online) PVB order correlation length f for 
different cylinder circumferences C at J2 = 0.3. Black circles 
indicate YC cylinders; red diamonds indicate XC cylinders. 
The green straight line illustrates the expected linear behavior 
at the quantum critical point of the two-dimensional system. 
The inset shows the local PVB order parameter vs. distance 
from the end of the cylinder in lattice spacings. The estimated 
PVB correlation length is indicated for each cylinder. 

XC9-1 and XC12-0. (We consider cylinders which don't 
accommodate PVB order, such as XC8-0 and XClO-0, in 
Appendix A.) We pin the PVB pattern at the cylinder 
ends by the choice of which spins are kept and how long 
the cylinder is (see Fig. 2]). 

We define a local PVB order parameter at each site 
using the nearest-neighbor spin correlations on the 3 ad- 
jacent plaquettes. The plaquettes form 3 sublattices in 
the PVB phase, as labeled in Fig. |3] (a)-(b), and one 
plaquette from each sublattice is adjacent to each site. 

Let Ea = —J^Ai'^i ' ^j) ^^'"^ O'^^'" the 6 bonds 

around plaquette A. Then we define the local PVB order 
parameter as 

27r 47r 
P = Ea + Eb exp {—i) + Ec exp i—i) . (2) 

Near the ends of the cylinders, this order parameter is 
nonzero, and in general it is a complex number. We 
extrapolate this local order parameter versus the trun- 
cation error to estimate its values in the limit of large 
bond dimension to. For long cylinders, its magnitude 
decays exponentially with distance from the end, with 
a correlation length that depends on the cylinder, as 
shown in the inset of Fig. [SJ These PVB correlations can 
be slightly incommensurate, particularly for the narrower 
YC cylinders. For the XC9-1 cylinder with a shifted con- 
nection, we measure the distance from the end along the 
direction perpendicular to the wrapping vector, to obtain 
the shortest PVB correlation length ^p. 

The PVB correlation length versus cylinder circum- 
ference C is shown in Fig. [5] This figure includes cylin- 
ders with all orientations, and we see that for our larger 
circumferences, the XC cylinders appear to have a longer 
PVB correlation length than the YC cylinders. If the 
2D system is at a quantum critical point, the correlation 
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FIG. 6: (color online) Spin singlet (solid line with filled sym- 
bols) and triplet (dashed line with open symbols) gaps for 
various cylinders versus inverse of circumference at J2 = 0.3. 

length is expected to be proportional to the circumfer- 
ence, by standard finite-size scaling. It appears that the 
correlation length actually increases faster than the cir- 
cumference, suggesting that this system may have weak 
PVB long range order in the 2D limit of infinite circum- 
ference. 

To explore the low-lying excited states on these cylin- 
ders, we apply DMRG to find the singlet and triplet gap. 
We lock the PVB pattern on the cylinder ends to guar- 
antee that the excited state doesn't include translations 
of the PVB pattern. For the triplet gap, we target the 
ground state in the 5*^ = and = 1 sectors separately. 
We make sure the states have the expected properties, 
e.g. the iSz = state has no signs of excitations and 
the 5*2 = 1 state has the excitation in the bulk, not at 
an edge. For the singlet gap, we target the ground state 
in the S'^ = sector first, until this state is very accu- 
rate. Then we target one more state in the same sector, 
with the sweeping restricted to exclude the end regions. 
The advantage is that the singlet excitations will always 
appear in the center of cylinder. We perform the calcu- 
lation on long cylinders to make sure that the gap is well 
converged. We use length 20 for XC cylinders and length 
30 for YC cylinders. 

The triplet and singlet gaps for XC and YC cylinders 
are shown in Fig. [S] The triplet gap decreases as the 
cylinders get wider, with the YC cylinders having higher 
triplet gaps than the XC cylinders. The triplet excita- 
tions on the XC cylinders appear to consist of two well- 
separated spinous. This effect is more pronounced when 
the cylinder gets wider. In YC cylinders, it appears that 
the spinous remain tightly bound. We don't yet under- 
stand why the behavior is so different between the XC 
and YC cylinders. From the trend of the triplet gap with 
circumference, particularly for the YC cylinders, we es- 
timate that in 2D limit the triplet gap is no more than 
0.12. The singlet gaps also show substantial variation 
between different cylinders. From the trends with cir- 
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cumference, we estimate that in the 2D hmit the singlet 
gap is no more than 0.07. 

In conclusion, we have studied the S — 1/2 honeycomb 
Ji — J2 Heisenberg model on various cylinders extensively 
using DMRG. We find that the ground state displays a 
paramagnetic phase for 0.26 ^ J2 ^ 0.36. By studying 
PVB order on various cylinders, we find that the PVB 
correlation length grows at least linearly with the cylin- 
der circumference. This suggests that in this phase the 
system is either quantum critical or has weak long-range 
PVB order. 
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FIG. 7: The ground state energy per site at J2 = 0.3 for var- 
ious cylinders with circumferences C. The two dashed lines 
are E/N — —0.4250 and —0.4260, which are the approximate 
upper and lower bounds of ground state energy that we esti- 
mate for the infinite two-dimensional system. 
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Appendix A: Ground state energy per site for 
J2 = 0.3. 



In this section, we present the ground state energies 
for J2 = 0.3 for various cylinders in Fig. [T) The energy 
per site for a given cylinder is calculated by subtracting 
the energies of two long cylinders with different lengths 
3l|. We find that the energy is much lower on the nar- 
row YC4-0 cylinder (C = 6), probably due to the con- 
tribution to the energy from short resonating paths that 
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FIG. 8: Top panel shows the ground state of a XC8-0 cyhnder 
for J2 — 0.3. It has a spinon with Sz — 1/2 and an anti-spinon 
with Sz = —1/2 on the edges. The bottom panel shows the 
excited state in the Sz ~ 1 sector with two spinons on the 
edges. These states have almost exactly the same energy. 

wind around the cylinder. This effect persists to a lesser 
extent for YC6-0 (C ~ 9). On all the other cylinders, 
the energy is consistently around -0.4255(5) with wider 
cylinders having larger uncertainty due to the extrapola- 
tion to large m. Thus we estimate an approximate upper 
bound on the energy as -0.4250. This bound is lower 
than other recent bounds from the variational entangled 
pair states method {E < -0.4210) 19] [? ] and VMC 
(i; < -0.4169) [13. 



Appendix B: Even and odd topological sectors on 
XC8-0 and XClO-0 cylinders. 

In this section, we will discuss the XC8-0 and XClO-0 
cylinders in detail. These cylinders are not compatible 
with PVB order, and show interesting behavior. The 
frustration of the PVB order leaves two different nearly- 
degenerate ground states corresponding to different topo- 
logical sectors. The sectors can be selected by adding 
or removing sites from the left and right edges. For 
zigzag-column left and right edges, the XC8-0 and XCIO- 
ground states have two free spinons on the cylinder 
edges, leading to four nearly degenerate states (Fig. [8]). 
The edge spinons can be thought of as unpaired spins 
in a resonating valence bond picture. Nearest-neighbor 
dimer coverings of a cylinder are in two distinct topolog- 
ical classes, specified by whether a vertical cut through 
the cylinder breaks an even or odd number of dimers. We 
can remove one site from each end of an XC cylinder to 
eliminate the end spinons. Since there are then N-1 sites 
remaining on the edge of XCN cylinder, this forces the 
dimer cover to be in the odd sector. 

The opposite topological sector also can be produced 
by adding two vacancies inside the cylinder. The effect 
of vacancies is to force the the opposite topological sec- 
tor between the two vacancies. The two vacancies need 
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FIG. 9: (color online) (a) Initial state for the even topological 
sector of the XC8-0 cylinder at J2 = 0.3. The width of the 
lines are proportional to \{Si ■ Sj) — e\ with e=0. A vertical 
line cuts through two singlet bonds, (b) The VBC ground 
state for the even sector with e=-0.32. (c) The spin-liquid- 
like ground state for the odd sector with e=-0.29 for horizontal 
bonds and e=-0.34 for vertical bonds. The spin correlations 
are fairly uniform in the center. On this XC8-0 cylinder the 
absolute ground state is that of the odd sector. 



to sit in special locations and not too far apart; other- 
wise spinons may bind to the vacancies, thus keeping the 
system in the same sector, (not shown) 

For the XC8-0 cylinder, the higher energy sector is the 
even topological sector, which is related to a VBC with 
an 8 site unit cell with a hexagon and a dimer, shown 
in the middle panel of Fig. [9] To get the energy per 
site for the even sector, we cut the ends and strengthen 
dimer bonds to favor the even sector initially (Fig. IHltop 
panel) . Then we remove these strengthened bonds to run 
the DMRG to get the final states (Fig. [H] middle panel). 
The energy splitting for the XC8-0 cylinder is shown on 
the left two panels of Fig. [TlJ The odd sector energy is 
always lower in the intermediate phase. This state looks 
more uniform (Fig. [Hlbottom panel) than the even sector, 
which shows strong valence bond order, but it has signif- 
icant anisotropy. In the intermediate phase regime, the 
energy splitting increases with J2, then decreases as J2 
approaches the second phase transition point. We didn't 
include results from higher J2 values, since the energy 
splitting is already small compared to the uncertainty in 
the extrapolation of the energy. The fact that the energy 
splitting reduces with J2 as one approaches the second 
phase transition point means that the phase transition in 
to the SVBC phase happens at essentially the same J2 
value in either sector. In Fig. [3] the odd sector was used, 
but we also did the same analysis in the even sector, ob- 
taining the same estimates of J2 at the phase transition. 
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FIG. 10: (color online) Similar to previous figure, (a) and (b) 
are the initial and the ground state for the even topological 
sector of the XClO-0 cylinder with J2 = 0.3. (c) is the odd 
topological sector. Unlike the XC8-0 cylinder, the even sector 
with apparent valence bond order is actually the ground state. 
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FIG. 12: (a) is XC8-0 cylinder with delta bond applied 
from 0.02 (left edge) to -0.05(right edge), with positive 
delta strengthen the vertical bonds and negative delta value 
strengthen the horizontal bonds. Thus we cut the cylinder 
edge to favor odd sector for positive delta and even sector for 
negative delta. Spinous locate between even and odd sectors, 
(b) is XClO-0 cylinder with delta from 0.01 to -0.01. 
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FIG. 13: Delta value versus sum of Sz along the cylinder 
zigzag direction for XC8-0 and XClO-0 cylinder. The data 
are extracted from previous figure. The left side of spinon 
{S < 0) is even sector, with odd sector on the spinon right 
side. We could see that for ground state with 5 = 0, XC8-0 
is in odd sector and XClO-0 is in even sector. 



FIG. 11: (color online) Different topological sector energy 
splittings for the XC8-0 and XClO-0 cylinders versus J2 with 
error bars, (a) and (c) are the energy splittings for the XC8-0 
and XClO-0 cylinders respectively, (b) and (d) are the energy 
per site for both even and odd sectors on XC8-0 and XClO-0 
cylinders. It's clear that on the XC8-0 cylinder the odd sector 
energy is lower. On the XClO-0 cylinder, the energy splitting 
is quite small, and becomes consistent with zero at J2 = 0.34. 



The even sectors of XClO-0 cylinder are shown on Fig. 
[TUl The even sector has a 20 sites VBC unit cell. The 
final state has stronger valence bond order than the XC8- 
cylinder. As opposed to XC8-0 cylinder, we show in the 
right two panels of Fig Hi] that the even sector is lower 
in energy. We could measure the spin correlation length 
from placing staggered magnetic fields on the edge and 
measuring how Sz decays. The odd sector has a much 



longer spin correlation length with ^ 12, as compared 
with roughly 3 in even sector. A longer spin correlation 
length usually indicate smaller spin triplet gap. We find 
that spin triplet gap is roughly 10"'^ for the odd sector 
compared with 0.075 for the even sector. Thus these 
results indicate that on XClO-0 cylinder the ground state 
is a valence bond crystal with short range spin correlation 
and a nonzero spin triplet gap. 

To further check that the even sector has lower energy, 
we strengthen and weaken some bonds to favor both sec- 
tors on one cylinder; one on the left, the other on the 
right, see Fig. [121 A domain wall with spinon will reside 
between these different sectors. Depending on whether 
the spinon is more pushed into the even or the odd sector 
at 6 = 0, we can determine which sector is the ground 
state. The results are consistent on both cylinders, shown 
in Fig. [13 
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Appendix C: Entropy analysis for intermediate 
phase with J2 — 0.3. 

In this section, we will discuss the von Neumann en- 
tanglement entropy (EE) in detail for J2 = 0.3 for all the 
XC and YC cylinders. For a cylinder with a non-shifted 
connection, eg. the XCN-0 and YCN-0 cylinders, we cal- 
culate the EE for splitting the lattice into subsystems A 
and B. The partition is a vertical line cut through the 
cylinder, that always breaks the same number of hori- 
zontal or diagonal bonds, i.e. breaking N/2 bonds for 
XCN-0 cylinder and N bonds for YCN-0 cylinder. We 
didn't consider the EE for cylinders with shifted connec- 
tions. 

In a 2D disordered gapped system, the EE obeys the 
area law [32] with 



S = aL + o{l/L) 



(3) 



where a is a non-universal constant, L is the length of the 
boundary between regions A and B. The second term dis- 
appears as L — 00. For a gapped topologically-ordered 
SL phase without broken symmetries, a extra term is in- 
cluded 



S = aL + "f + o{l/L) 



(4) 



with 7 is the topological EE[33|, [SJj. For a minimum 
entangled state with Z2 topological order, j — — ln(2) at 
T = 0. Thus the EE is a constant independent of length 
for a fixed width cylinder, since the boundary length is 
always the same. For a 2D critical gapless system, the 
EE scales as: 



S = aL + b\nL + c{L) ln[sin( — )] + o(l/L) (5) 

L 



for a, L X L torus [35|-|37| 



The second term comes from 
X is subsystem size and the 



gapless Goldstone mode 
third term is the function depends on "chord length". 
Thus the EE depends on system width L and subsys- 
tem size X for a 2D gapless system, but independent of 
subsystem size for a gapped system. 

In Fig. [TiT a). we show the EE versus subsystem size 
X for YCN-0 cylinders. Since YC cylinders have strong 
plaquettes appearing on the cylinder edges, we would ex- 
pect that the EE should oscillate with a PVB pattern, 
with two high values and one low value corresponding to 
cuts at strong and weak bonds of plaquettes. With the 
decay of PVB order to the center, we also observe the 
decay of the EE oscillation pattern and measure the EE 
when its value saturates to a constant value. We notice 
that YC6-0 has a longer PVB correlation length than 
YC4-0, since the EE oscillates with more periods near 
the YC6-0 cylinder edges. On the YC8-0 cylinder, the 
state in which we measure the EE is not well converged, 
and thus it has stronger PVB order. The EE shows an 
apparent PVB order oscillation pattern. 
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FIG. 14: (color online) The EE versus subsystem size x for 
(a) YC cylinders, (b) XC6-0 and XC12-0 cylinders, (c) even 
and odd sector in XC8-0 cylinders, (d) even and odd sec- 
tor in XClO-0 cylinders. Note on YC8-0 cylinder, the EE 
is measured for a state with m = 5000. This state is not 
well converged and it has larger PVB order as seen in the EE 
oscillation pattern. 
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FIG. 15: Saturated EE versus cylinder circumference C for 
the XC and YC cylinders. 



For the XC cylinders in Fig. [T4lb-d). the entropy is 
saturated in the center, independent of subsystem size for 
the ground state. For the XC8-0 cylinder, the odd sector 
EE is independent of the sub-system size x and has the 
same value for cylinders of different lengths. The EE has 
much stronger dependence on the sub-system size for the 
even sector. On the contrary for the XClO-0 cylinder, 
the EE saturates in the center of the cylinder for the 
even sector. All these results again show that the GS 
for all the cylinders with J2 = 0.3 are gapped with short 
spin correlation lengths. 

In Fig. [151 we present the EE versus cylinder width 
to extrapolate the topological EE 7. It is surprising that 
not all the data points collapse onto a line. The EE 
from XC6-0 and XClO-0 deviates from a linear extrap- 
olation. It is possible that XC6-0 cylinder is too nar- 
row, so that J2 = 0.3 is too close to the second phase 
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transition point. The XClO-0 cylinder is a valence bond 
crystal with 20 sites unit cell in the even sector. We find 
that 7 = —0.619, which is close to gapped Z-i SL with 
7 = In 2 = -0.693. 

In a recent paper by Hong-Chen Jiang et al. [lOj , the au- 
thors extrapolate the EE to get the TEE with 7 ~ — In 2 
and identify Z2 SL phase on many models. As those au- 
thors realize, for the above method to be accurate when 
identifying topological phases, all the correlation lengths 
(spin-spin, dimer-dimer, etc.) have to be short compared 



with the cylinder widths, so that DMRG results will al- 
ways lead to a minimum entangled state (MES) with zero 
vison flux through the cylinder. For the model in this 
paper, even though the spin-spin correlation length is 
short, PVB correlation lengths get longer as the cylin- 
der gets wider. Thus the extrapolation method for the 
TEE would lead to strong finite size effects. Therefore 
the non-zero 7 we find in this paper doesn't indicate that 
the GS for the Heisenberg J\ — J2 model on the honey- 
comb lattice is a Z2 SL in the 2D limit. 



